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Abstract 


A  thorough  investigation  of  leading  edge  heat  transfer  on  a  model 
geometry  has  been  performed  using  Large-eddy  simulation  (LES),  with 
support  from  the  current  and  a  previous  AFOSR  grant.  The  results  from 
the  leading  edge  study  are  presented  as  a  detailed  Appendix  in  this  report. 
To  address  the  problem  of  solving  the  full  three-dimensional  turbulent 
flow  in  a  turbine  passage,  including  the  effects  of  free-stream  turbulence 
from  the  combustor,  and  the  unsteady  interactions  between  various  blade 
rows  in  the  form  of  moving  wakes  and  shocks  it  was  necessary  to  extend 
the  LES  solver.  A  high-order  overset  LES  code  suitable  for  the  study  of 
turbomachinery  passages  has  been  developed.  This  approach  is 
summarized  in  this  report.  Some  preliminary  results  have  been  obtained 
using  this  code,  mostly  as  cases  for  code  validation.  The  code  is  currently 
being  used  to  investigate  flow  in  turbine  passges  under  realistic  conditions, 
and  detailed  comparisons  with  available  experiments  will  be  made. 
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1  Introduction  and  Objectives 


The  current  trend  in  design  of  gas  turbine  engines  is  to  achieve  a  reduction  in  cost  and  weight 
of  the  engine  and  at  the  same  time  not  compromise  on  the  thrust  developed  by  them.  There  are 
primarily  two  avenues  to  achieve  this  viz (i)  by  using  fewer  blades  and  heavily  loading  them*  and 
(it)  by  increasing  the  turbine  inlet  temperature.  As  a  result  the  nozzle  guide  vanes  and  rotor  blades 
of  a  high-pressure  turbine  (HPT)  are  exposed  to  a  hostile  thermal  environment.  Active  cooling  of 
the  blade  surfaces  becomes  necessary  in  order  to  maintain  their  integrity.  This  has  a  performance 
penalty  since  the  air  required  to  cool  the  blades  is  bled  from  the  compressor.  The  increased  loading 
of  the  blades  leads  to  boundary  layer  seperation  and  results  in  performance  losses.  Therefore  heat 
transfer  prediction  of  the  thermally  critical  regions  and  flow  prediction  in  the  blade  passage  are  a 
vital  elements  of  the  overall  improvement  of  gas  turbine  engines. 

The  primary  objective  of  this  project  was  to  investigate  fundamental  issues  controlling  aerodynamic 
losses  and  elevated  heat- transfer  to  turbine  vane  and  rotor  blades  using  large  eddy  simulations 
(LES).  LES  simulation  methodology  was  used  to  study  the  physics  behind  simplified  problems 
relevant  to  turbomachinery. 

(a)  Simulation  of  realistic  free-stream  turbulence  with  intensity  and  length  scale  appropriate  for 
turbine  blade  heat- transfer  and  study  the  heat- transfer  augmentation  in  idealized  geometry  and  its 
interaction  with  blade  leading  edge  region 

(b)  Simulation  of  the  disturbed  turbulent  boundary  layers  over  a  realistic  turbine  blade  surface  {in 
a  blade  row  cascade)  in  the  presence  of  frec-strcam  turbulence  and  combustor  hot  streaks  and  the 
resulting  blade  heat  transfer  distribution 

Other  flow  processes  such  as  blade  row  interactions,  tip  vortices,  tip  clearance  flows,  and  secondary 
flows  in  the  blade  row  that  are  critical  to  the  overall  predictions  in  a  realistic  turbine  blade  row 
were  not  considered. 

A  high  order  overset  LES  code  was  developed  in  order  to  perform  numerical  investigation  of  the 
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above  listed  problems.  In  the  following  sections,  we  will  look  at  the  form  of  the  Navier-Stokes 
equations  suitable  for  the  problems  of  interest  and  discuss  the  implementation  in  the  LES  code. 
We  will  then  look  at  some  results  obtained  using  this  high  order  LES  code. 


2  Governing  Equations 

The  equations  governing  the  flow  of  a  compressible  Newtonian  fluid  obeying  Fourier's  law  of  heat 


conduction  can  be  written  in  terms  of  the  contravariant  velocity  components  (lA),  in  a  general 
curvilinear  coordinate  system  (x1)  as  (Axis  (1983)) 

~  +  (pvkU  =  0  0) 

I \(PV')  +  =  -  [(p  +  s'J]  +  (sr’V*  +  s'Vfe)]  J  ( 2 ) 

^.+  [(B+p)vi]j  -  +  {^9ikvk).  (3) 

where,  gtj  and  g'3  are  the  covariant  and  contravariant  metric  tensors  respectively.  The  equations 
are  recast  in  the  following  form  to  facilitate  numerical  implementation 

+  WkMpvt)  "  0  <4) 

+  +  =  g|i  {v^  9IJ(p+5/.»fi)  } 

+  y/g  |pij'  (p  +  -  ( 5 ) 


j^E)  +  ±[^E  +  p^\  =  |j{v4S«|;+rW]}  (6) 

where  J  =  ^/g  is  the  Jacobian.  Hence,  a  natural  choice  of  the  dependent  vector  variables  would  be 

q  =  { y/gp ,  y/gPv\JgE}, 

A  filter  <7(x)  with  compact  support  is  then  applied  to  the  governing  equations  in  the  domain  of 
interest  (fi)  to  derive  the  equations  for  LES.  Any  tensor  quantity  A1^ '  ^  may  be  filtered  as 

Vf4±(x,«)  =  JJJnAi£±{x',t)G{x-x')y/gdxl 
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The  curvilinear  equivalent  of  Favre  filtering  can  be  written  as 


itl—tn 


..  .  JgA\  , 

J<1— >m  _ 

jl"in  y/9P 

Using  this  definition,  we  can  write  the  filtered  governing  equations  as 


|w&) +  £<*&•*>  -  o 


dxJ 

d»^ij  p 


a 

9iJ 


r« 


ftp*  *  Q)  *9S 

a  a  a 


5j(\/5  £)  +  gjj  [\/s(£  +  p)  ®*]  =  - fcj  (^)  +  aS  »  8‘) 

~9iJ  ( 

where  the  resolved  and  sub-grid  scale  stress  tensor  and  heat  flux  vector  are 

fV  -  + 

/iePr  <9x* 

sg.s'T^  =  \/PP  (v‘W  -  i>V) 

*9^J  =  (TV  -  TV) 


In  addition,  the  equation  of  state  becomes 


VSP  = 


7  -  1 


\f9PT 


(?) 


(8) 


(9) 


(10) 

(11) 

(12) 

(13) 


(14) 


The  resulting  sub- grid  scale  (SGS)  stress  tensors  are  closed  using  the  dynamic  Smagor insky  model 
(See  Nagarajan  (2004)  for  the  details). 


3  Numerical  Implementation 


The  flow  in  a  turbine  passage  can  be  split  into  two  regions,  namely  (i)  a  near  wall  region  where 
the  turbulent  boundary  layer  lias  to  be  resolved,  (ii)  and  the  region  further  away  from  the  blade 
surface  where  primarily  the  evolution  of  the  free-stream  turbulence  has  to  be  raptured.  This  implies 
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a  wry  fine  grid  requirement  in  the  boundary  layer,  while  the  region  farther  away  imposes  a  far 
less  stringent  grid  requirement.  To  simulate  this  kind  of  flow,  an  overset  grid  approach  is  suitable, 
where  the  near  wall  region  uses  a  O-type  bodyfitting  grid  which  is  embedded  in  a  periodic  H-type 
grid  in  the  passage.  Such  a  capability  has  been  implemented  in  the  code  written  by  Nagarajan  et 
al  (2003), 


3.1  Discretization  and  time  marching  schemes 


A  sixth-order  compact  finite  difference  scheme  of  Lele  (1992)  is  used  to  discretize  the  variables  in 
space*  A  staggered  arrangement  of  variables  is  employed,  where  the  scalar  variables  are  stored  at 
the  center  of  the  grid  celts  while  the  velocity  components  are  stored  at  their  respective  cell  faces. 


In  order  to  optimise  the  cost  of  the  computations,  a  mixed  implicit-explicit  time  marching  algorithm 
is  used*  The  bodyfitting  O-grid  region  is  marhed  in  time  using  a  Beam- Warming  type  implicit 
scheme 


3£/*+1  -4Un  +  Un-{ 
2A  t 


=  -G(C/n+1,*n4-*) 


(15) 


where  the  solution  at  time  levels  n  and  (n  —  1)  are  used  to  find  the  solution  at  time  level  (n  +  1). 
In  order  to  reduce  the  cost,  of  the  inversion,  an  approximate  factorization  of  the  right  hand  side 
Gn+l  is  applied.  The  background  H-type  mesh  is  marched  using  the  standard  RK-3  scheme.  See 
Nagarajan  (2004)  for  the  details  of  the  implementation. 


3.2  Intergrid  Communication 

The  soul t ions  in  the  two  grid  zones  are  coupled  through  an  overlap  zone,  where  the  data  is  in¬ 
terpolated  between  the  grids  using  the  fourth  order  Her  mite  interpolation  scheme  developed  by 
Dclfs  (2001).  In  two  dimensions,  the  Hermit.e  interpolation  at  an  overlap  point  is  constructed  us¬ 
ing  the  function  values  and  the  first  derivatives  at  the  four  surrounding  points  which  bound  the 
interpolated  point  in  the  uniform  computational  space*  A  local  grid  system  (£,f?)  is  constructed  in 
the  computational  space,  with  the  origin  of  the  coordinate  system  at  the  center  of  the  box  formed 
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by  the  points  bounding  the  interpolated  point.  The  coordinates  of  the  interpolated  point  in  this 
coordinate  system  are  found  by  the  inverse  mapping  (£,7?)  ==  The  forward  mapping 

is  defined  at  all  points  using  the  Hermite  interpolation  scheme,  and  the  inverse  mapping 
A/_1{z,  y)  is  found  using  a  Newton-Rapbson  procedure.  The  interpolation  formula  is  give  by 


m,v)=  Y,ci^)ft+u+i'  +  cfk(z,T,) 


Idf 


t,k=0 


Lao 


i+(,j+fc 


+  C|fc(Ci  *?) 


df 


Lar, 


i+lj+k 


+  0(  A4)  (16) 


where  Cjj^,  Cfk,  and  are  coefficients  of  interpolation.  The  details  of  the  implementation  can 
be  found  in  Delfs  (2001).  Whereas  Delfs  uses  the  Dispersion  Relation  Preserving  (DRP)  finite 
difference  scheme  for  computing  derivatives,  we  use  the  sixth  order  compact  finite  difference  scheme. 


4  Summary  of  LES  Code  development 

A  high-order  parallel  compressible  LES  code  developed  by  A,  Xiong  (2003)  was  used  to  study 
leading  edge  heat  transfer.  Another  high-order  parallel  compressible  LES  capability  to  simulate 
transitory  flows  was  developed  under  DoD  ASC  support  by  S,  Nagarajan,  based  on  the  algorithm  by 
Nagarajan,  Lele  and  Ferziger  (2003).  Nagarajan  et  al  (2007)  used  this  solver  to  study  the  problem 
of  bypass  transition  on  model  geometries  with  focus  on  the  effects  of  leading-edge  bluntness.  It 
was  possible  to  adapt  any  of  these  codes  to  do  the  full  three-dimensional  turbomachinery  passage 
simulations.  Nagarajan’s  code  was  chosen  for  the  robustness  and  efficiency  of  his  algorithm  and  its 
implementation. 

Under  this  grant,  Nagarajan’s  LES  code  was  extended  to  study  the  flow  and  heat  transfer  in  full 
turbine  cascade  passages.  As  a  first  step,  iso- thermal  wall  boundary  condition  was  implemented  to 
study  heat  transfer.  The  conflicting  requirements  of  maintaining  turbulence  resolution  in  the  tree- 
stream,  capturing  the  vortex  dynamics  of  turbulence  as  it  is  distorted  in  the  turbine  passage  and 
especially  around  the  leading-edge  region,  maintaining  sufficient  resolution  in  the  viscous  boundary 
layers,  and  allowing  the  substantial  flow  turning  representative  of  turbine  cascades  required  the 
development  of  a  overset  LES  capability. 


5 


A  computational  issue  in  the  simulation  of  HP  turbine  flows  is  the  treatment  of  moving  shocks 
which  may  arise  in  these  transonic  flow  fields.  Recently*  Fiorina  &  Lele  (2007)  have  developed 
a  high-order  nonlinear  hyper-viscosity/hyper-diffusivity  based  shock  capturing  scheme  for  LES  of 
supersonic  combustion  flows*  such  as  a  jet  in  supersonic  cross-flow  in  SCRAM  jet  applications. 
These  flows  involve  complex  three-dimensional  shock  waves  and  its  interaction  with  turbulence. 
We  intend  to  implement  this  shock-capturing  methodology  in  the  overset  LES  code  in  the  near 
future. 

The  overset  mesh  algorithm  is  designed  to  run  efficiently  on  parallel  computers  and  for  this  purpose 
special  parallel  solvers  for  large- banded  matrices  were  implemented. 


5  Results 

The  effect  of  free-stream  turblence  on  stagnation  point  flow  and  heat  transfer  was  studied  using 
the  LES  capability  developed  by  Xiong  (2003).  It  was  observed  that  in  the  presence  of  free-stream 
turbulence*  intense*  quasi-streamwise  vortices  develop  near  the  leading  edge  region*  which  were 
responsible  for  heat  transfer  augmentation.  The  results  from  this  study  are  currently  under  review 
for  publication  in  JFM,  The  manuscript  has  been  included  as  an  appendix  to  this  report* 

The  new  overset  LES  code  is  being  used  for  the  turbine  passage  simulation.  This  overset  LES 
capability  has  been  tested  on  2D  and  3D  cases.  The  DNS  of  flow  past  a  circular  cylinder  at 
Re.  =  150  and  Ma  =  0*2  was  simulated  on  a  2D  mesh.  A  body-fitted  O-type  grid  of  size  225  x  170 
and  a  background  H-t.ype  grid  of  size  670  x  500  was  used*  Figure  1  shows  the  dipole  nature  of  the 
sound  radiated  by  the  cylinder.  Figure  2  shows  the  excellent  comparison  with  the  pressure*  lift  and 
drag  coefficients  obtained  from  this  simulation  and  the  2D  DNS  by  Inoue  and  Hatakeyama  (2002). 
A  typical  two  dimensional  simulation  with  half  a  million  grid  points  took  aboiit  12  hrs  on  8  Intel 
Xeon  processors  in  order  to  get  converged  results. 

Exploratory  2D  and  3D  simulations  using  the  C3X  HP  vane  geometry  were  performed.  Figures  3-4 
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show  results  from  a  2D  DNS  of  the  C3X  HP  cascade*  A  combination  of  O-type  grid  of  size  400  x  45 
around  the  blade  and  a  background  H-type  grid  of  size  400  x  175  in  the  passage  was  used*  The 
inlet  Reynolds  number  based  on  the  inlet  velocity  and  chord  of  the  blade  was  50,000  and  the  inlet 
Mach  number  was  0.15.  The  simulation  was  initialized  by  the  solution  from  a  RANS  computation 
of  the  passage  using  an  in-house  RANS  code. 


A  preliminary  3D  LES  calculation  was  performed  using  the  same  grid  as  the  2D  calculation,  but 
with  the  2D  planar  grid  stacked  in  the  spanwise  direction  to  form  a  3D  grid.  32  spanwise  points 
were  used  over  a  spanwise  length  of  30%  chord*  The  solution  from  the  2D  calculation  was  used 
to  initialize  the  3D  calculation.  Although  no  free  stream  turbulence  was  imposed,  the  wake  from 
the  blade  was  however  observed  to  become  turbulent*  A  way  to  quantify  aerodynamic  losses  in 
the  passage  is  to  the  compute  the  wake  loss  coefficient  along  a  line  in  the  wake.  The  wake  loss 
coefficient  is  defined  as 


Q  = 


Pt>in  Pj,cx 
Pt,in  -  Ps,  ex 


(17) 


The  wake  loss  profile  for  the  3D  LES  of  the  C3X  cascade  is  shown  in  figure  5.  The  values  seem  to 


compare  well  with  the  experiments  of  Ames(1994)]  although  the  Reynolds  number  of  the  calculation 
was  much  lower  and  there  was  no  external  free  stream  turbulence. 


These  exploratory  calculations  have  allowed  us  to  gain  confidence  with  the  new  Overset  LES  code. 
Simulations  matching  experimental  conditions  and  with  imposed  extremal  free-streani  and  wake 
turbulence  will  be  conducted  in  the  near  future* 


6  Summary 

A  thorough  investigation  of  leading  edge  heat  transfer  has  been  performed  using  LES,  with  support 
from  the  current  and  a  previous  AFOSR  grant*  The  results  from  the  leading  edge  study  have  been 
presented  at  conferences  and  have  also  appeared  as  articles  in  archival  journals  (see  references 
section).  The  focus  has  then  shifted  towards  solving  the  full  three-dimensional  problem  of  flow 
in  a  turbine  passage,  including  the  effects  of  freostream  turbulence  from  the  combustor,  and  the 
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unsteady  interactions  between  various  blade  rows  in  the  form  of  moving  wakes  and  shocks,  A  high- 
order  overset  LES  code  suitable  for  the  study  of  turbomachinery  passages  has  been  developed.  Some 
preliminary  results  have  been  obtained  using  this  code,  mostly  as  cases  for  code  validation.  The 
code  is  currently  being  used  to  investigate  flow  in  turbine  passges  under  realistic  conditions,  and 
detailed  comparisons  with  available  experiments  will  be  made.  This  study  is  expected  to  continue 
under  a  future  grant  from  AFOSR. 

7  Personnel 

Professor  Sanjiva  Lele  is  the  principal  investigator  for  this  project.  A  former  graduate  student  Dr, 
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Figure  1;  (a)  The  grid  used  for  the  cylinder  simulation,  showing  the  overlap  region  between  the 
O-  and  H-type  grids,  (b)  Fluctuation  pressure  field  AP,  showing  the  dipolar  nature  of  sound  field 
radiated  by  the  cylinder.  The  contour  levels  are  from  APmin  =  -0.1  M25  to  APmai  =  0.1  M2 5 
with  an  increment  of  0,0025M2  5.  solid  line,  AP  >  0;  -  ■  AP  <  0. 
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Figure  2:  Cylinder  simulation:  (a)  Time  averaged  pressure  coefficient,  (b)  Time  history  of  lift 
and  drag  coefficients.  Lines:  Current  simulation.  Symbols:  Simulation  by  I  none  and  Hatakeyama 
(2002). 


Figure  3:  2D  C3X  HP  cascade  simulation:  (a)  The  grid  used  for  the  C3X  HP  cascade  simulation, 
showing  the  overlap  region  between  the  O-  and  H-type  grids.  Every  fourth  grid  point  has  been 
shown,  {b}  pv  contours  showing  the  wake  shed  by  the  blade. 
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Figure  4:  2D  C3X  HP  cascade  simulation:  (a)  Pressure  and  (b)  heat  transfer  to  the  blade. 


Figure  5:  3D  C3X  HP  cascade  simulation:  (a)  Iso-surfaces  of  spanwise  vorticity  show  the  turbulent 
nature  of  the  wake,  (b)  Wake  loss  coefficient  along  a  line  in  the  wake  perpendicular  to  the  blade 
chord. 
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In  this  paper,  the  effects  of  free-stream  turbulence  on  stagnation  point  flow  and  heat 
transfer  are  investigated  through  Large  Eddy  Simulation  (LES)  of  homogeneous  isotropic 
turbulence  impinging  upon  an  isothermal  elliptical  leading  edge.  Turbulent  mean  flow  and 
Reynolds  stress  profiles  along  the  stagnation  streamline,  where  the  mean  flow  is  strain 
dominant,  and  at  different  downstream  locations,  where  the  mean  flow  gradually  becomes 
shear-dominated,  are  used  to  characterize  evolution  of  the  free-stream  turbulence.  The 
Reynolds  stress  budgets  are  also  obtained,  and  the  turbulence  anisotropy  is  analyzed 
through  the  balance  between  the  mean  flow  strain  and  the  velocity  pressure  gradient 
correlation.  In  the  presence  of  free-stream  turbulence,  intense,  quasi -stream wise  vortices 
develop  near  the  leading  edge  with  a  typical  diameter  on  the  order  of  the  local  boundary 
layer  thickness.  These  strong  vortices  cause  the  thermal  fluxes  to  peak  at  a  location 
much  closer  to  the  wall  than  that  of  the  Reynolds  stresses,  resulting  a  greater  sensitivity 
to  free-stream  turbulence  for  the  heat  transfer  than  the  momentum  transfer.  The  heat 
transfer  enhancement  obtained  by  the  present  LES  agrees  quantitatively  with  available 
experimental  measurements.  The  present  LES  results  are  also  used  to  examine  the  eddy 
viscosity  and  pressure  strain  correlations  in  Reynolds  stress  turbulence  models. 
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1.  Introduction 
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Stagnation  point  flow,  where  fluid  approaching  a  solid  surface  divides  into  diverging 
streams,  occurs  ubiquitously  in  nature  as  well  as  in  many  engineering  problems.  An 
improved  understanding  of  the  stagnating  flow  turbulence  and  heat  transfer  is  critical  to 
a  wide  range  of  engineering  applications,  from  protecting  a  gas  turbine  blade  from  being 
melted  by  the  extremely  hot  and  turbulent  combustion  gas,  to  enhancing  the  efficiency 
of  modern  micro- electronics  cooling  system  using  micro-jet  array  impingement. 

In  his  pioneering  work,  Hiemenz  established  that  stagnation  point  flow  is  one  of  the 
very  few  types  that  admit  exact  solutions  of  Navier-Stokes  equations.  Subsequent  studies 
have  been  focused  on  understanding  perturbed  stagnation  point  flows,  i.e.  the  evolution 
and  effects  of  various  disturbances,  including  turbulence,  that  are  present  in  the  free 
stream.  One  of  the  prominent  effects,  first  observed  in  1920s  {Piercy  k  Richardson  1928, 
1930),  is  that  free-stream  turbulence  induces  large  heat  transfer  enhancement  at  the 
stagnation  point  (Giedt  1949;  Hegge-Zijnen  1957;  Kestin  et  a l.  1961).  Various  semi- 
empirical  correlations  have  been  proposed  to  quantify  this  effect,  mostly  based  on  the 
mean  flow  Reynolds  number  and  the  free-stream  turbulence  intensity  {Smith  k  Kuethe 
1966;  Kestin  k  Wood  1971;  Lowery  k  Vachon  1975;  Mehendale  et  at .  1991),  but  they 
are  met  with  limited  success.  Later,  turbulence  length  scale  was  recognized  to  play  a 
critical  role  in  determining  the  overall  effects  of  free-stream  turbulence.  By  incorporating 
turbulence  length  scale  into  the  correlations,  the  applicability  of  prediction  correlation 
is  significantly  improved  (Ames  k  Moffat  1990;  Van  Fossen  et  a l.  1995;  Dullenkopf  k 
Mayle  1995).  Although  different  definitions  of  length  scale  are  used  in  these  correlations, 
the  heat  transfer  enhancement  is  found  generally  to  increase  with  increasing  Reynolds 
number  and  turbulence  intensity,  but  decrease  with  increasing  turbulence  length  scale. 

Theoretical  studies  of  the  free-stream  turbulence  effects  in  stagnation  point  flows  have 
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ranged  from  the  mean  flow  unsteadiness,  the  linear  /nonlinear  instability,  to  the  vortex 

st  retch  ing /am  p  I  i  fi  c  ation  effects.  Light  hi  II  (1954)  studied  a  pulsating  mean  flow  around  a 
cylindrical  body  and  obtained  the  Stokes- layer  correction,  but  no  significant  change  in 
heat  transfer  was  found.  Similarly  modulated  stagnation  point  flows  were  also  studied 
by  Fshigaki  (1970),  Pedley  (1972)  and  Merchant  k  Davis  (1989),  but  with  an  emphasis 
on  the  skin  friction.  The  linear  stability  analysis  was  initiated  by  Gortler  (1955)  and 
Hammerlin  (1955),  further  extended  by  Kestin  k  Wood  (1970),  and  finally  clarified  by 
Wilson  k  Glad  well  (1978).  It  has  shown  that  plane  stagnation  point  flow  is  always  lin¬ 
early  stable  to  three  dimensional  disturbances.  For  finite  amplitude  disturbances,  Lyell  k 
Huerre  (1985)  showed  that  the  flow  can  be  destabilized  if  the  level  of  the  external  two-  or 
three  dimensional  disturbances  exceeds  a  threshold  value.  Recently,  the  linear  instability 
for  the  more  general  attachment  dine  bound  ary  layer  flow  has  also  been  investigated  by 
Lin  k  Malik  (1996)  and  Theofilis  et  ai  (2003).  Morkoviu  (1979)  argued  that  the  enhance¬ 
ment  of  heat  transfer  is  more  likely  a  result  of  forced  response  to  upstream  disturbances 
as  opposed  to  any  internal  flow  instability,  a  view  advocated  earlier  by  Sutera  (1965)  who 
analyzed  the  amplification  of  incoming  organized  disturbances  by  mean  strain,  indicat¬ 
ing  the  sensitivity  of  the  heat  transfer  to  vortical  disturbances.  For  turbulence  impinging 
on  a  circular  cylinder  with  either  very  large  or  very  small  scales,  Hunt  (1973)  obtained 
the  behavior  of  second  order  turbulence  statistics  by  using  a  generalized  rapid  distor¬ 
tion  theory  (RDT).  Xiong  k  Lele  (2004)  showed  the  critical  parameter  in  determining 
the  evolution  of  upstream  three-dimensional  disturbance  is  the  ratio  between  the  distur¬ 
bance  length  scale  to  the  Hiemenz  boundary  layer  thickness.  Using  numerical  simulation 
of  a  swept  Hiemenz  boundary  layer,  S  pal  art  (1989)  found  that,  out  of  initial  white- noise 
disturbances,  the  most  unstable  disturbance- mode  is  the  one  with  the  same  similarity 
form  as  the  mean  Hiemenz  flow,  an  assumption  made  in  the  stability  analysis  mentioned 
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above.  The  flow  structures  induced  by  free-stream  turbulence  in  a  stagnation  region  were 
found  to  be  qualitatively  similar  to  those  induced  by  upstream  organized  disturbances 
(Xiong  &  Lele  2001).  Bae  et  at,  (2000)  showed  that  different  length  scales  generate  quite 
different  flow  patterns  and  in  turn  different  heat  transfer  responses  in  a  plane  stagnation 
point  flow.  Although  much  progress  has  been  made  over  the  years,  a  complete  under¬ 
standing  of  the  effects  of  free-stream  disturbances,  particularly  turbulence,  in  stagnation 
point  flows  has  not  been  achieved. 

The  free- stream  turbulence  effects  have  also  proven  difficult  to  be  incorporated  in  en¬ 
gineering  turbulence  models.  The  one  equation  model,  e.g,  Spalart  &  Allmaras  (1992), 
while  widely  used  and  shown  to  be  particularly  successful  in  aerodynamic  flows  {Bardina 
et  aL  1997;  Wilcox  2001),  does  not  explicitly  account  for  the  effects  of  free-stream  turbu¬ 
lence.  Standard  two  equation  models,  e.g,  k  -  t  or  k  —  ui  model,  when  used  in  stagnation 
point  turbulent  flows,  badly  over  predict  the  turbulent  kinetic  energy  and  beat  trans¬ 
fer  —  this  is  termed  as  “stagnation  point  anomaly”  (Durbin  1996;  Champion  U  Libby 
1994,  1991).  A  fundamental  difficulty  with  the  two  equation  models  is  that  turbulence 
is  assumed  to  be  predominately  isotropic  and  in  quasi-equilibrium.  But  stagnating  flow 
turbulence  can  be  strongly  anisotropic  because  different  components  of  the  fluctuation 
velocity  respond  differently  to  the  mean  flow  straining.  Reynolds  stress  models  provide 
the  generality  to  account  for  the  anisotropy  of  the  turbulence  explicitly.  Im  et  ol.  (2002) 
used  three  variants  of  Reynolds  stress  model,  the  GL  model  (Gibson  Launder  1978}, 
the  GL-GL  model  (Craft  et  at.  1993)  and  the  SSG  model  (Speziale  et  aL  1991),  to  com¬ 
pute  both  the  impinging  and  counter- current  stagnation  flows.  Although  the  predictions 
From  these  models  were  better  than  A;  —  €  type  models,  all  the  models  still  severely  over¬ 
predict  the  turbulence  kinetic  energy,  and  have  large  discrepancies  in  other  Reynolds 
stress  components  when  compared  to  experimental  measurements.  The  problem  stems 
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from  the  over  prediction  of  the  energy  production  and  the  under  prediction  of  the  redis¬ 
tribution  by  the  pressure  strain  correlations.  This  uncertainty  may  be  attributed  to  the 
difficulty  in,  and  hence  the  lack  ofT  detailed  experimental  measurements  in  the  immediate 
vicinity  of  the  leading  edge. 

In  this  study,  we  carry  out  LES  of  the  impingement  of  free- stream  turbulence  upon  an 
isothermal  elliptical  leading  edge.  To  allow  direct  comparison  with  experimental  results, 
the  flow  configuration  and  simulation  parameters  are  taken  to  match  the  experiments 
by  Van  Fossen  et  a L  (1995).  The  first  goal  of  the  present  study  is  to  gain  an  improved 
understanding  on  the  evolution  of  stagnating  flow  turbulence  and  the  mechanism  of  the 
heat  transfer  enhancement.  Secondly,  the  simulation  is  aimed  to  help  reduce  some  of  the 
uncertainty  in  turbulence  modelling  for  strain-dominated  flows.  The  relative  magnitudes 
and  distributions  of  various  turbulence  quantities  are  obtained  from  the  LES  data.  The 
present  simulation  shows  that  the  largest  change  in  the  turbulence  structure  occurs  in 
the  immediate  region  of  the  leading  edge.  Thirdly,  it  is  hoped  the  present  LES  data 
may  help  the  general  development  and  calibration  of  turbulence  models  when  turbulence 
anisotropy  becomes  strong. 

The  rest  of  the  paper  is  organized  as  follows.  The  basic  governing  equations  for  the  LES 
are  given  in  section  2.  The  numerical  methods  and  validation  procedures  are  described  in 
section  3,  along  with  a  blending  scheme  for  generating  free-stream  turbulence.  The  main 
simulation  results,  including  the  mean  flow  profiles,  turbulence  intensity  and  Reynolds 
stress  budget,  are  presented  in  section  4,  In  section  5,  we  summarize  the  results  and  give 
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Figure  1 .  Flow  configuration  and  inflow,  outflow  and  wall  boundary 

2.  Navler-Stokes  equation  and  SGS  modeling 

The  flow  configuration  for  the  present  LES  is  shown  in  figure  i  with  both  the  laboratory 
coordinates  {x,y}z)  and  the  local  surface  coordinate  labelled.  The  mean  flow  is 

two-dimensional  and  the  turbulence  is  assumed  to  be  homogeneous  in  spanwise  direction 
z * 

The  governing  equations  for  the  LES  of  compressible  turbulence  are  the  filtered  Navier- 
Stokes  equations.  With  the  dynamic  subgrid  stress  {SGS)  model  for  compressible  flows 
{Moin  et  aL  1991),  they  are  written  as 

P.t +  (?«*).<  =  o  (2.1) 

=  ~P,i  +  ^  +  C2AT-5*i)TjJ  (2.2) 

pf,t + puif,i  +  (7  -  lypfujj  =  + 

7(7  ~Rl)M  [iSiiSji  +  2jiSi}§ij],  (2.3) 

where 

P=^3  +  IC,?Ai|S|!  (2'4) 


At  =  A  +  CpA2]S|fle 


(2,5) 
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kT  =  k  +  CpA2\S\^^-.  (2.6) 

7  Fn 

Here  the  tilde  stands  for  Farve  average;  p  is  the  density;  u,  is  the  velocity  vector,  and 
T  is  the  temperature.  pT  is  the  subgrid  eddy  viscosity,  kt  is  the  subgrid  eddy  thermal 
conductivity,  and  Prt  is  the  turbulent  Prandtl  number.  S,j  =  +  uJtl)  is  the  rate-of- 

strain  tensor.  M ,  Re  are  the  mean  flow  Mach  number  and  Reynolds  number,  and  in  the 
present  study  they  are  based  on  the  upstream  mean  velocity  t/^,  the  speed  of  sound 
and  the  leading  edge  diameter  D.  Pr  =  0.7  is  molecular  Prandtl  number  and  7  =  1.4  is 
the  ratio  of  specific  heats.  The  eddy  coefficients  C\  for  the  subgrid  normal  stress  ts  set  to 
zero  and  the  turbulence  Prandtl  number  Prt  is  set  to  unity.  The  eddy  coefficient  C  for 
the  subgrid  shear  stress  is  computed  using  the  standard  dynamic  SGS  model,  see  (Xiong 
2004}  for  details. 

3,  Numerical  Method 

The  numerical  method  and  validation  problems  for  the  present  LES  are  described 
briefly  in  this  section.  For  more  details,  the  reader  is  referred  to  Xiong  (2004), 

3X  Implicit  scheme  with  linearized  subiterations 
We  first  recast  the  compressible  flow  governing  equations  in  a  general  form: 

Ut  +  F(U)  =  0  (3.1) 

where  U  =  {pT  u7  vt  w,  T}T  is  the  vector  of  flow  variables  and  F(U)  represents  the  nonlin- 
ear  and  viscous  terms.  Since  for  wall-bounded  flows,  explicit  time  integration  schemes  are 
prohibitively  expensive  because  of  the  CFL  stability  constraint,  in  this  study  we  use  an 
implicit,  dual-time  stepping  scheme  with  linearized  subiteration,  which  may  be  expressed 
as 

7  A  ^ 

[/  +  -f  Ar.4(£/n)]A(/k  =  -Ar7Jk 


(3.2) 


8 

with 


Z.  Xiong  and  S.  K.  Lele 


A 


dF 

dV' 


nk  =  wk  4^un~l  +  mk)> 


a  uk  =  -  uk . 


(3  3) 


Here  At  is  the  physical  time  step,  and  Ar  is  the  pseudo  time  step  for  subiteration.  At 
each  physical  time  step,  Un  is  taken  to  be  the  initial  value  for  Uk  to  start  the  subiteration 
at  k  —  0.  If  the  subiteration  converges,  i.e,  A Uk  — *  0  one  gets  £/fc+1  =  Uk .  The  final 
value  of  £/fc+1  is  taken  to  be  £/n+1,  and  ftrt+1  “  0  recovers  the  second  order  fully  implicit 
scheme,  which  is  unconditionally  stable.  Notice  that  instead  of  a  function  of  Uk  as  in 
standard  subiteration  schemes,  the  left  hand  side  operator  in  (3.2)  is  only  a  function 
of  £/n;  therefore  (3.2)  is  linear  for  variable  A A  LU  decomposition  of  the  coefficient 
matrices  is  performed  at  the  first  step  of  each  subitereation,  and  the  factored  matrices 
are  stored  and  used  until  the  subiteration  converges.  This  obviates  the  need  for  inverting 
the  coefficient  matrices  at  every  subiteration  step,  and  therefore  significantly  improves 
the  subiteration  efficiency. 


3.2.  Spatial  Discretization 

At  interior  nodes,  the  five-point  central  difference  scheme  is  used  for  first  and  second 
derivatives.  Near  the  computational  boundaries,  five- point  biased  stencils  are  used.  The 
resulting  difference  schemes  are  fourth  and  third  order  accurate  for  the  first  and  second 
derivatives,  respectively.  Owing  to  the  intrinsic  non-dissipative  nature,  central  difference 
schemes  are  often  subject  to  the  so-called  two-  6  wave  instability  caused  by  the  decoupling 
of  the  even  and  odd  grid  points.  In  the  LES  context,  one  should  not  expect  that  a  physics- 
based  SGS  model  will  automatically  suppress  such  an  instability.  Beside  the  fact  that 
the  SGS  model  itself  is  often  implemented  with  the  same  central  difference  scheme,  the 
length  scale  of  the  two- 5  waves  is  much  shorter  than  that  of  the  SGS  model.  The  purpose 
of  the  SGS  model  is  to  extract  energy  at  a  correct  rate  at  a  cutoff  scale —  the  scale  at 
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which  the  discrete  representation  becomes  inadequate —  and,  for  a  grid  spacing  hy  this 

is  typically  around  a  wavenumber  of  hk  w/2.  But  for  the  two-(5  wave,  it  occurs  at  the 

Nyquist  wavenumber  hk  —  7 r.  Hence  to  suppress  the  two-5  waves,  appropriate  numerical 

damping  procedures  must  be  applied.  In  the  present  simulations,  we  choose  the  following 

fourth-order  dissi pat i on 


(3.4) 


where  cj  is  the  amplitude  of  the  dissipation,  and  A€,  A,,  and  A*  are  the  grid  spacing  in 
streamwise,  normal  and  spanwise  directions  in  computational  space.  Care  are  must  be 
taken  to  ensure  that  the  added  numerical  dissipation  is  minimal  and  does  not  deteriorate 
the  resolved  solution.  For  this  purpose,  e <t  is  chosen  to  be  such  that  the  magnitude  of 
Dn  is  significantly  smaller  than  the  truncation  error  of  the  difference  schemes.  Based  on 
modified  wave  number  analysis,  this  may  be  achieved  by  requiring  ^  O.OIA^,  where 
N3  =  (N^  +  Nn  +  Nz)/ 3  is  the  average  number  of  grid  points  in  one  spatial  direction.  The 
resulted  Dn  is  one  order  of  magnitude  lower  than  the  truncation  error  and  its  effects  on 
the  results  of  the  present  LES  is  expected  to  be  negligible.  A  more  detailed  discussion  of 
the  method  and  the  demonstration  of  its  effectiveness  can  be  found  in  Xiong  (2004). 


3.3.  Boundary  Conditions 

Boundary  conditions  are  introduced  to  replace  the  governing  equations  at  the  inflow, 
outflow  and  the  wall  boundaries  of  the  computational  domain.  Consider  an  arbitrary 
boundary  constraint  at  time  level  n  +  1  on  the  flow  variable  U  —  {pfuf  u, 


B(Un+*)  =  0. 


(3.5) 
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The  genera]  implicit  treatment  of  the  boundary  conditions,  in  terms  of  the  variable 
6U  —  Un+l  —  t/n,  can  be  written  as 


(3  6) 


The  boundary  condition  at  the  inflow  provides  the  upstream  mean  flow  information  and 
also  introduces  free-stream  disturbances  or  turbulence  into  the  computational  domain. 
For  a  subsonic  inflow,  characteristic  analysis  shows  that  four  incoming  quantities  must  be 
specified  along  with  one  outgoing  quantity  computed  from  the  interior  domain.  The  spe¬ 
cific  choice  of  these  quantities  depends  on  the  formulation  of  the  problem.  In  this  study, 
we  constrain  the  entropy,  spanwise  and  tangential  velocities,  and  the  incoming  Riemann 
invariant.  The  outgoing  Riemann  invariant  is  computed  by  first  order  extrapolation  from 
the  interior  points  close  to  the  boundary. 

At  the  outflow,  the  parabolized  Navier-Stokes  equations  is  used  following  Collis  (1997), 
i.e.  the  stream  wise  second  order  derivatives  in  the  equations  are  neglected.  In  addition, 
the  pressure  gradient  at  the  outflow  is  obtained  from  the  corresponding  potential  flow 
solution.  This  treatment  has  shown  to  yield  an  adequate  and  stable  outflow  boundary 
conditions  both  for  the  laminar  and  turbulent  flow  computations.  No  slip  and  isothermal 
conditions  are  applied  at  the  wall  for  the  velocities  and  temperature,  and  periodicity 
condition  is  imposed  in  spanwise  direction. 

Prior  to  turbulence  simulations,  the  code  is  first  validated  on  problems  of  laminar 
compressible  boundary  layer  at  a  leading  edge.  The  velocity  and  temperature  profiles  are 
compared  with  analytic  self-similar  solutions.  An  leading  edge  acoustic  receptivity  prob¬ 
lem  is  also  computed  and  compared  with  previous  numerical  studies.  These  validation 
results  are  summarized  in  Appendix  A  ♦ 
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3.4.  Generation  of  free- stream  turbulence 

In  this  section,  we  describe  a  blending  procedure  combining  independent  but  statistically 
equivalent  realizations  of  homogeneous  isotropic  turbulence  into  a  unified  turbulent  flow 
field  which  serves  as  a  realistic  representation  of  free-stream  turbulence.  Consider  two 
such  turbulence  fields*  f/^3^  and  that  need  to  be  concatenated  in  x  direction,  we 

introduce  a  blending  zone,  extending  from  the  interface  into  each  of  the  fields  by  a  short 
distance,  and  a  linear  combination 

u  =  cos#  u*1*  -f  sin#  u*2\  (3.7) 

where  #  varies  smoothly  from  0  to  tt/2  across  the  blending  zone.  The  new  field  varies 
smoothly  from  one  field  to  another  and  retains  the  mean  values  and  second  order  statistics 
of  the  original  fields.  The  dependence  of  #  on  x  within  the  blending  zone  introduces 
an  extra  dilatation  field,  but  it  can  be  removed  by  using  the  Helmholtz  decomposition 
theorem  for  the  velocity  vector.  In  figure  2,  we  plot  the  energy  spectra  of  two  original 
and  the  blended  turbulence  fields.  It  can  be  seen  that  the  energy  spectrum  of  the  blended 
field  remains  essentially  the  same  as  the  originals.  A  a  very  small  amount  of  energy  is 
present  at  the  lowest  wave  number  due  to  the  increased  length  in  the  x-direction.  More 
comparisons  of  energy  spectra,  as  well  as  other  quantities  of  interest,  can  be  found  in 
Xiong  et  al.  (2004), 

In  the  present  LES,  twelve  realizations  of  such  independent  but  statistically  identical 
homogeneous  isotropic  turbulence  fields  are  prc-computed  using  LES  with  the  same  dy¬ 
namic  SGS  model.  To  keep  the  same  spanwise  domain  size  as  in  the  main  LES,  these 
simulations  are  carried  out  in  a  rectangular  box  of  the  size  1.6D,  1  .CD  and  0,4Z?  in  x,  y 
and  z  direction,  respectively.  A  128  x  128  x  32  grid  is  chosen  to  ensure  the  isotropy  of 
the  resulting  turbulence.  The  code  is  adapted  from  a  DNS  code  {Lui  2003)  which  uses 
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FlCUHE  2.  Energy  spectra  for  the  original  and  blended  data  sets. - :  data  1, 

- data  2, - blended 

a  sixth-order  compact  finite  differencing  scheme  and  a  fourth  order  Runge-Kutta  time 
integration.  Once  the  twelve  turbulence  fields  are  obtained j  they  are  lined  up  spatially 
and  joined  together  at  the  interfaces  by  applying  the  above  blending  procedure.  The 
resulting  turbulence  field,  twelve  times  longer  than  each  individual  realization  but  with 
the  same  characteristics*  serves  as  the  free-streara  turbulence,  which  will  be  convected 
into  the  computational  domain  through  the  inflow  boundary. 

4,  LES  results 

4.1.  Simulation  conditions 

The  flow  condition  for  the  present  LES  is  taken  from  the  wind  tunnel  experiments  by 
Van  Fossen  et  al  (1995),  and  specifically,  corresponding  to  the  data  set  No.  244.  In  the 
experiments,  free- stream  turbulence  is  generated  by  placing  a  turbulence- generating  grid 
upstream  of  a  leading  edge  model.  The  shape  of  the  leading  edge  is  a  half-ellipse  and 
the  aspect  ratio  (length  ratio  between  the  major  and  minor  axis)  is  3  for  this  data  set. 
Different  values  of  turbulence  intensity  and  length  scale  are  obtained  by  varying  the 
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Reo 

Tu 

L/D 

Rgl 

Rct 

Tw/T0 

domain  (z,y,  z) 

N,  x  Nv  x 

42,000 

0.15 

0.05 

0.1 

210 

93 

1.075 

(3.5D,  5D,  0.4D) 

191  x  144  x  48 

Table  1.  Parameters  of  the  mean  flow  and  free-stream  turbulence  for  the  present  LES, 


mesh  size  of  the  grid  and  its  location  relative  to  the  leading  edge.  The  model  is  kept 
at  a  constant  temperature,  and  the  heat  transfer  rate  is  measured  in  the  stagnation 
point  region.  The  mean  flow  and  free- stream  turbulence  parameters  are  summarized  in 
Table  1.  The  mean  flow  Reynolds  number  He#  —  U^Djv  is  based  on  far  upstream 
incoming  velocity  Uoo  and  the  leading  edge  diameter  of  curvature  D.  Moo  —  UoofQ<x>  is 
the  mean  flow  Mach  number  based  on  the  sound  speed  far  upstream.  Tu  and  L  are 
the  free-stream  turbulence  intensity  and  integral  length  scale;  they  are  determined  at  the 
location  of  the  leading  edge  but  in  the  absence  of  the  model.  The  numerical  procedure 
for  determining  L  from  the  auto-correlation  measurements  is  given  in  Van  Fossen  et  al . 
{1995}  and  applied  in  the  present  LES.  Rcl  —  TuL/v  is  the  Reynolds  number  based 
on  integral  length  scale  and  turbulence  intensity  and  the  Rtr  is  the  Taylor  microscale 
Reynolds  number,  Tm  and  T0  are  the  surface  temperature  on  the  wall  and  the  total  or 
stagnation  gas  temperature  in  the  free-stream. 

The  computational  grid  is  generated  by  an  algebraic  multi-surface  method  (Eiseman 
1985)  in  the  x  —  y  plane  which  guarantees  the  grid  orthogonality  at  the  wall  and  the 
inflow /outflow  boundaries.  The  grid  points  are  clustered  towards  the  wall  and  the  stag¬ 
nation  point  but  uniformly  spaced  in  spanwise  direction.  In  the  present  LES  study,  a 
dual-grid  approach  is  used.  First,  a  large  outer  grid  is  generated  covering  half  of  the  el¬ 
lipse  downstream  and  extending  vertically  from  the  bottom  to  the  top  wind  tunnel  wall. 
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A  potential  flow  that  matches  the  experimental  base  flow  conditions  is  computed  on  this 
grid.  For  the  subsequent  LES  calculations*  a  smaller  inner  grid  covering  only  the  leading 
edge  region  is  extracted  from  the  outer  grid.  On  this  inner  grid,  the  the  inflow /out  flow 
boundary  conditions  for  the  mean  flow  can  therefore  be  provided  by  the  potential  flow 
solution.  This  dual- grid  approach  enables  us  to  improve  the  grid  resolution  near  the  lead¬ 
ing  edge  at  a  reduced  computational  cost.  For  the  LES  grid,  the  minimum  stream  wise 
grid  spacing  is  A xmin/D  —  CL 0037  at  the  stagnation  point,  and  the  minimum  normal 
spacing  is  AymiTt/Z?  =  0.00058  at  the  wall  The  Hiemenz  scale  for  the  laminar  boundary 
layer*  y/v/A/D  —  0.0036*  is  resolved  by  7  grid  points  in  the  normal  direction,  here  A  is 
the  outer  flow  strain  rate.  The  time  step  for  the  present  LES  is  At  —  0.0015,  and  the 
total  integration  time  AT  At  for  obtaining  the  turbulence  statistics  is  25. 

4.2.  Mean  flow  results 

The  contour  plots  for  the  mean  flow  streamwise  velocity  U,  transverse  velocity  V  and 
temperature  T  are  shown  in  figure  3.  Outside  the  boundary  layer*  the  difference  be¬ 
tween  the  mean  and  the  laminar  velocity  contours  are  minimal*  and  along  the  stagnation 
streamline,  the  velocity  profiles  are  almost  indistinguishable  (not  shown  here*  see  Xiong 
(2004)).  It  confirms  the  experimental  measurements  that  the  stagnation- line  velocity  is 
essentially  independent  of  the  free-stream  turbulence.  This  is*  however*  not  the  case  for 
the  temperature  profile.  As  shown  in  figure  4a,  a  steeper  gradient  is  developed  at  the 
wall  for  the  mean  temperature  profile  (about  25  percent  higher  than  the  laminar  value)* 
indicating  the  heat  transfer  enhancement  in  the  presence  of  free-stream  turbulence. 

The  comparisons  between  the  mean  and  laminar  profiles  of  the  tangential  velocity 
Us  and  the  normalized  temperature  (T  -  1  )/(Tw  -  I)  are  further  shown  in  figure  4b-d* 
using  the  local  $  -  n  coordinates  at  three  downstream  locations  s/D  =  0.2*  0.8  and  1.6 
(marked  in  figure  3c).  The  change  of  the  slope  at  the  wall  is  again  about  25  percent  for  the 
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Figure  3,  Turbulent  mean  velocity  and  temperature  contours.  Re  =  42,000,  Mo  —  0,15* 
Tu  —  0.05,  L/D  —  0.1,  a)  streamwise  velocity  U,  contour  minimum:  0.0  ,  maximum:  1.3* 
increment:  0.05,  b)  transverse  velocity  V,  minimum:  -0.6,  maximum:  0.6t  increment:  0,05.  c) 
temperature  T,  minimum  0.995,  maximum  1.075,  increment:!). 005.  The  streamwise  locations  of 
the  four  cross  sections  are  also  marked. 
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Figure  4.  Comparison  between  turbulent  mean  and  laminar  velocity  and  temperature  profiles 

at  different  streamwise  locations, - :  turbulent  mean,  — - :  laminar.  s/D  =  0, 02, 0,8, 1.6 

from  a)  to  d). 

temperature  and  only  about  5  percent  for  the  velocity,  signaling  a  greater  sensitivity  to 
the  free-stream  turbulence  for  the  heat  transfer  than  for  the  skin  friction.  Moreover,  the 
profiles  of  the  mean  streamwise  velocity  also  show  that,  despite  the  impinging  turbulence, 
the  boundary  layer  remains  pre- transitional,  A  longer  streamwise  distance  is  needed  than 
the  present  simulation  domain  before  a  fully  turbulent  boundary  layer  can  develop. 


4,3.  Vortex  dynamics 

Figure  5  shows  a  visualization  of  instantaneous  spanwise  velocity  within  the  x  -  y  plane 
at  z  =  0,  and  the  temperature  gradient  on  the  model  surface  near  the  stagnation  point 
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FIGURE  5,  Instantaneous  span  wise  velocity  in  x  —  y  plane  and  temperature  gradient  on  the 

body  surface. 

s  =  0.  It  can  be  seen  that  the  isotropic  free-stream  turbulent  eddies  have  been  strongly 
stretched  in  the  stream  wise  direction  near  the  stagnation  point.  The  ensuing  impinge¬ 
ment  of  these  vortices  on  the  leading  edge  modifies  the  underlying  thermal  boundary  layer 
dramatically.  On  the  model  surface ,  the  resulting  heat  transfer  distribution  develops  into 
thin,  streamwise  elongated  streaky  structures.  The  corresponding  instantaneous  temper- 
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Figure  6-  Instantaneous  temperature  contours  at  different  locations  on  the  body  surface, 

s/D  -  —0,8  -  0.8* 

ature  fields  are  shown  in  figure  6  through  a  series  of  wall  normal  sections  at  different 
streamwise  locations.  The  temperature  contours  show  distinctive  mushroom-like  struc¬ 
tures.  By  tracing  these  structures  in  consecutive  sections,  one  easily  sees  the  elongated 
streaky  vortices  being  wrapped  around  the  leading  edge  by  the  mean  flow. 

To  better  understand  the  typical  vortex  structures  and  their  effects  on  heat  transfer, 
the  temperature  contour  and  the  corresponding  velocity  field  in  the  stagnation  plane 
s  “  0  are  shown  in  figure  7,  The  velocity  fields  show  clearly  that  strong,  amplified 
y-oriented  vortices  produce  reverse  flow  (u  <  0)  in  the  stagnation  region,  which  lifts  up 
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Figure  7.  Temperature  contour  and  the  corresponding  velocity  field  in  the  stagnation  plane 
at  different  times.  The  mean  flow  direction  is  downward  and  the  time  interval  is  A t  -  0.6. 


the  hot  fluid  dose  to  the  wall  and  swaps  it  with  the  lower  temperature  fluid  away  from 
the  wall  In  doing  so,  the  velocity  disturbance  creates  the  mushroom-like  structures  in 
the  temperature  contours.  Directly  underneath  these  mush  room -like  structures,  however, 
the  thermal  boundary  layer  is  thicker  than  the  undisturbed  ease  —  the  heat  transfer  rate 
decreases.  It  is  in  the  region  between  these  mushroomdike  structures,  where  the  disturbed 
flow  has  a  normal  velocity  directed  towards  the  wall,  that  the  boundary  layer  become 
thinner,  and  consequently  the  heat  transfer  increases.  The  overall  spanwise  averaged  heat 
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transfer  rate  is  thus  determined  by  the  distribution  and  intensity  of  these  thickened  and 
thinned  regions. 

As  shown  in  figure  7,  the  mushroonvtike  structures  typically  do  not  remain  at  a  fixed 
spanwise  location;  instead  they  move  laterally  over  a  significant  distance  during  their 
lifetime.  This  indicates  the  lateral  movement  of  the  strong  streamwise  vortices  by  which 
these  mushroom  structures  are  formed.  When  a  vortex  approaches  a  wall,  vorticity  is 
generated  at  the  wall  through  the  no-slip  boundary  condition;  therefore,  the  lateral  move¬ 
ment  of  the  vortices  may  be  characterized  through  the  movement  of  the  wall  vorticity. 
By  following  the  slowest  descent  line  on  the  space-time  correlation,  the  average  lateral 
speed  for  the  wall  vorticity  along  the  stagnation  line  (z-axis)  is  found  to  be  ±0.06,  which 
is  about  half  of  the  maximum  wrms  at  the  leading  edge.  The  mechanism  of  the  lateral 
movement  can  be  understood  through  the  so-called  wall-blocking  effect.  Considering  a 
single  vortex  approaching  a  wall,  an  inviscid  image  vortex  with  opposite  sign  is  induced 
at  an  equal  distance  on  the  other  side  of  the  wall  to  enforce  the  no-penetration  wall 
boundary  condition.  The  mutual  induction  of  these  two  vortices  generates  a  tangential 
velocity  whose  magnitude  is  proportional  to  the  strength  of  the  vortices.  This  causes  the 
lateral  movement  of  the  approaching  vortex  in  the  spanwise  direction,  and  in  turn  the 
movement  of  the  mushroom  structures  in  the  temperature  contours. 

The  stretching  of  the  free-stream  turbulence  eddies  may  further  be  characterized  quan¬ 
titatively  using  spanwise  two-point  correlations  of  velocity.  Figure  8  shows  the  streamwise 
velocity  correlation  at  different  locations  along  the  stagnation  streamline.  At  the  inflow 
boundary,  the  correlation  length  is  long,  representing  the  length  scale  of  the  free-stream 
turbulence  in  the  absence  of  the  strain  effect.  Once  the  turbulence  enters  the  domain,  the 
spanwise  length  scale  decreases  monotonically  and  dose  to  the  leading  edge  it  reaches 
a  scale  of  local  boundary  layer  thickness.  A  model  problem  on  the  competition  between 
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Figure  8,  Two  point  correlation  of  streamwise  velocity  in  spanwise  direction  along  stagnation 

streamline. . :  x  -  1.55D, - :  x  —  2.0 £>, - — :  x  =  2.97 D  (in  the  free  stream),  and 

- :  s  =  2.99D  (within  the  boundary  layer)* 

the  amplification  due  to  vortex  stretching  and  damping  due  to  viscosity  in  a  disturbed 
Hiemenz  boundary  layer  was  studied  by  Xiong  &  Lele  (2004 )T  and  the  leading  edge  region 
of  the  present  flow  represents  the  same  process  in  the  context  of  a  complete  flow  field.  We 
should  also  remark  that  because  of  this  mono  tonic  decrease  of  the  turbulence  spanwise 
length  scale,  the  choice  of  the  spanwise  simulation  domain  size  is  a  balance  between  en¬ 
compassing  the  large  scale  ffee-stream  turbulence  at  the  inflow  and  sufficiently  resolving 
the  intense,  small-scale  streamwise  vortices  at  the  leading  edge.  In  the  present  LES,  we 
place  our  focus  at  the  leading  edge  region  and  have  chosen  the  spanwise  domain  size  to 
be  four  times  the  integral  length  of  the  free-stream  turbulence.  A  similar  LES  with  span- 
wise  domain  size  about  seven  times  of  the  turbulence  length  scale  was  also  performed 
in  Xiong  (2004).  The  comparison  between  the  two  cases  suggests  that  the  LES  results 
are  not  sensitive  to  the  spanwise  domain  size  as  long  as  the  small  scale  vortices  at  the 
leading  edge  are  adequately  resolved. 

The  spanwise  averaged  non-dimensional  heat  transfer  coefficient,  denoted  by  Frossting 
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FrGURE  9,  Averaged  Frossling  number  distribution  on  the  body  surface.  - — :  turbulent 

mean, - ;  laminar,  o  :  experimental  data  (Van  Fossen  ef  ai.  1995}* 

number  Fr  =  Nu/^Reo ,  where  Nu  being  the  Nusselt  number,  is  plotted  in  figure  9 
as  a  function  of  the  streamwise  distance  from  the  stagnation  point.  The  experimental 
measurements  from  Van  Fossen  ei  at,  (1995)  are  also  plotted  For  comparison.  Both  the 
LES  and  experimental  results  show  that  the  shape  of  the  Fr  distribution  is  largely  the 
same  as  in  the  laminar  case,  but  their  values  are  significantly  higher  in  the  presence 
of  free-stream  turbulence.  Although  the  present  LES  result  is  slightly  lower  than  the 
experimental  measurements,  as  not  all  of  the  turbulence  scales  are  resolved  in  LES,  the 
overall  agreement  is  good.  As  shown  by  Xiong  &  Lele  (2004),  there  is  an  optimal  length 
scale  for  the  stagnation  point  heat  transfer  enhancement;  scales  larger  or  smaller  than  the 
optimum  produce  a  smaller  heat  transfer  enhancement.  In  order  to  accurately  predict  the 
free-stream  turbulence  effects,  it  is  thus  important  to  resolve  the  turbulence  structures 
up  to  the  optimal  length  scale.  For  this  reason,  the  effect  of  SGS  modelling  needs  to  be 
evaluated.  A  typical  distribution  of  the  instantaneous  SGS  eddy  viscosity  vT  is  shown 
in  figure  10a,  and  the  corresponding  time  averaged  distribution  is  shown  in  figure  10b. 
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FIGURE  10.  Contours  of  SGS  eddy  viscosity  vrjv.  a)  Instantaneous  (averaged  in  2),  contour  min¬ 
imum  :  0.05,  maximum:  0.75,  increment:  0.014.  b)  Turbulent  mean,  minimum  :  0.05,  maximum: 
0.52,  increment:  0.040. 

Both  figures  show  that  the  maximum  vr  occurs  in  the  leading  edge  region,  where  small 
scale  turbulence  is  produced  by  strong  vortex  stretching  and  amplification.  However, 
the  value  of  ur  is  relatively  small,  about  half  of  the  molecular  viscosity,  indicating  that 
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the  grid  resolution  close  to  the  wall  is  adequate  for  capturing  the  near  wall  eddies.  In 


this  sense,  the  present  LES  is  essentially  a  quasi- DNS  of  the  near  wall  flow.  However* 
because  of  the  grid  stretching  in  the  wall  normal  direction,  outside  the  boundary  layer 
the  turbulence  scales  are  not  completely  captured,  and  we  have  found  the  SGS  eddy 
viscosity  is  necessary  to  prevent  the  solution  from  diverging  caused  by  the  intermittent 
passage  of  strongly  stretched  vortices  outside  the  boundary  layer, 

4,4.  Reynolds  stress  and  turbulence  budgets 

Because  of  the  spanwise  homogeneity,  the  significant  Reynolds  stresses  in  the  present 
flow  are  ua,  v,2t  and  uV,  and  the  significant  heat  fluxes  are  uT',  v*V.  Along  the 
stagnation  streamline,  the  symmetry  in  y  direction  also  leads  to  urvf  =  vtTf  =  0, 

It  is  instructive  to  write  the  governing  equations  for  Reynolds  stress  and  heat  transport* 
Using  the  notion  of  Reynolds  average  and  letting  /  denote  the  time  averaged  value  of  /T 
we  may  decompose  the  turbulence  field  into  a  mean  and  a  fluctuation  part*  Le, 

p  =  p  +  p\  u,  =  Ui  +  txji  P  —  R  4-  p\  T  =  T  +  T\  (4,1) 
and  similarly  for  the  transport  coefficients 

A  =  A  +  A\  pr  —  p  +  *t  =  «  +  «'.  (4.2) 

The  steady  transport  equation  for  the  Reynolds  stress  tITuJ  can  be  written  as 

-  Vi.k  [p(uiUj-)  +  UkP%  +  ] 

“  Vj.k  (pK«!)  +  ukpfu\  +  ] 

+  (4-3) 

On  the  right  hand  side  of  (4.3),  the  first  row  is  the  velocity  pressure  gradient  correlation; 
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the  second  and  the  third  rows  are  the  production  terms;  the  fourth  row  is  the  transport  of 


Reynolds  stress  due  to  turbulent  fluctuation.  For  incompressible  flow,  it  can  be  reduced 
to  a  divergence  form  and  is  termed  as  turbulence  transport.  Here,  for  compressible  flow, 
we  use  the  term  turbulence  convection  to  refer  to  the  sum  of  these  two  terms.  The 
in  the  last  row  is  the  total  dissipation 


4>r  =  ^  [  (A  £>>'  +  D,i  X%  +  X'D'^r )  + 

(XD^  +  Dj  +  XDW)] 

(P5j*,*ui  +  Sjk,klt'K  +/i'5'iJku')] 

=  e  +  $„  (4.4) 


where  D  =  u1t  is  the  divergence  of  the  velocity  field;  e  is  the  turbulence  dissipation  rate 
defined  by 


E  — 


-Re^^ 


(4.5) 


and  the  remaining  term  $v  =  -  £  is  loosely  termed  as  viscous  diffusion  by  analogy 

to  the  incompressible  flow  case,  although  it  also  contains  terms  involving  dilatation  and 
fluctuation  in  the  turbulent  viscosity  coefficient. 

Similarly,  for  turbulence  heat  flux  u'7\  we  may  also  write 


pUkluTU  = 

-  T.jfc  \pu'ku\  +  Ukptu'i  +  p'u'ku\  ] 

-  Ui,k  [gr  +  ukW  + 


-  KKn*  -  VkA<T').k  -  Mwru 

+  p~^  + 


+  $p  +  $H 


(4.6) 


26  Z.  Xiong  and  $.  K.  Lele 


Figure  11.  Turbulence  intensity  along  the  stagnation  streamline,  — — :  urm.M . ■: 

On  the  right  hand  side  of  (4.6),  the  first  row  is  the  temperature  pressure  gradient  corre¬ 
lation;  the  second  and  the  third  rows  are  the  heat  flux  production;  the  fourth  row  is  the 
turbulence  convection  and  the  fifth  row  is  the  turbulent  heat  conduction.  In  the  last  row, 
$p  and  4>f/  are  the  pressure  dilatation  work  and  viscous  dissipation.  The  expressions  for 
4>p  and  4*//  are  omitted  here;  they  are  not  important  in  the  present  flow. 

4.5.  Turbulence  intensity  along  the  stagnation  streamline 
The  change  of  the  different  components  of  turbulence  intensity  characterizes  the  anisotropy 
of  stagnating  turbulence.  The  root- mean- square  value  of  velocity  fluctuation  along  the 
stagnation  streamline  is  shown  in  figure  11.  The  combination  of  the  the  turbulence  length 
scale  and  the  close  location  of  the  inflow  boundary  largely  eliminates  the  free  decaying 
stage  of  the  free-stream  turbulence.  After  a  relatively  balanced,  isotropic  development 
before  x  ~  2.4,  the  turbulence  develops  strong  anisotropy  between  x  ~  2.4  and  2.7, 
Within  this  region,  both  urms  and  wrTn,  increase  but  vrTn3  decreases,  which  is  consistent 
with  the  prediction  of  RDT  (Hunt  1973),  indicating  that  the  dominant  mechanism  for 
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FIGURE  12.  The  evolution  of  the  Reynold s- stress  invariants  along  the  stagnation  streamline, 
the  turbulence  anisotropy  is  by  vortex  stretching  (Lee  &  Reynolds  1985)  in  y  direction. 
At  closer  distances  x  >  2.7,  the  presence  of  the  model  surface  generates  a  pronounced 
inviscid  blocking  effect.  By  the  non- penetration  boundary  condition,  the  normal  velocity, 
Urmat  must  vanish  at  the  wall.  This  causes  urTn3  to  decay  over  a  distance  on  the  order  of 
the  turbulence  integral  length  scale,  typically  larger  than  the  boundary  layer  thickness. 
As  a  result,  the  kinetic  energy  of  urtn *  is  partially  redistributed  into  v  and  w  compo¬ 
nents,  leading  to  not  only  a  further  increase  in  tnrm31  but  also  an  increase  in  vrms  that 
was  previously  decreasing.  Within  an  extremely  close  distance  to  the  wall,  viscous  effect 
dominates  and  all  the  turbulence  vanishes  on  the  wall.  The  evolution  of  the  turbulence 
anisotropy  may  also  be  characterized  by  the  trajectory  in  the  plane  of  the  invariants  £ 
and  t 7,  here  £  -  (///*/ 2)1/3>  77  =  {-Z4/3}1/2,  and  Ilb  and  IIIb  are  the  second  and  third 
invariants  of  the  Reynolds  stress  anisotropy  tensor  6,j  =  tZJuj/u*Ufc  —  StJ/Z  (Pope  2000). 
This  is  shown  in  figure  12  by  the  discrete  dots  representing  different  locations  along 
the  stagnation  streamline.  It  can  be  seen  that  from  a  nearly  isotropic  state  the  Reynolds 
stress  anisotropy  develops  as  the  stagnation  point  is  approached  and  reaches  a  peak  value 
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7]  =  0.225  inside  the  boundary  layer.  The  turning  point  around  77  =*  0.1  corresponds  to 
the  peak  location  of  urms  in  figure  11,  demarcating  the  beginning  of  significant  inviscid 
blocking  effect. 

4.6,  Turbulence  budgets  along  the  stagnation  streamline 

The  budgets  for  the  normal  Reynolds  stresses,  u'2,  v12  and  wf2y  along  the  stagnation 
streamline  are  shown  in  figure  13.  In  these  figures,  the  total  dissipation  is  used  as  the 
dissipation  term.  For  un ,  the  production  term  -ut28Ujdx  is  associated  with  the  mean 
flow  strain  rate,  representing  an  energy  transfer  from  the  mean  flow  to  uf 2,  Conversely, 
the  production  for  va  is  v'2dV/dy ,  which  has  an  opposite  sign  and  represents  an  energy 
transfer  back  to  the  mean  flow.  There  is  no  production  term  for  wa  because  the  mean 
flow  is  two  dimensional.  The  turbulence  transport  has  the  largest  magnitude  in  ua  but 
is  confined  mostly  in  the  near  wall  region.  It  changes  sign  as  the  wall  is  approached, 
and  if  integrated  along  the  stagnation  line,  it  results  a  vanishing  net  contribution  to  the 
turbulence  energy.  Turbulent  convection  is  not  important  for  t/2  as  compared  to  other 
terms.  For  wt2t  we  found  that  the  turbulent  convection  term  largely  cancels  with  the 
mean  convection  term,  and  the  sum  of  these  two  exerts  a  weak  dissipative  effect  on  w'2. 
The  total  dissipation  extends  farthest  from  the  wall  for  uf2  and  closest  for  vf2.  It  has  the 
largest  value  for  u/2,  yet  the  smallest  for  u'2.  Except  for  the  dissipation  term  does 
not  play  a  significant  role  along  the  stagnation  line  except  extremely  dose  to  the  wall  for 
vf2. 

Of  particular  interest  for  the  stagnating  turbulence  is  the  redistribution  term,  i.e.  the 
correlation  between  velocity  and  pressure  gradient.  Along  the  stagnation  streamline,  it 
is  comparable  in  magnitude  for  all  three  normal  Reynolds  stresses.  In  strain- dominated 
flows,  its  primary  effect  is  the  generation  of  turbulence  anisotropy.  Being  largely  negative 
in  figure  13a  for  uf2,  it  is  positive  for  vf2  and  u/2,  as  shown  in  figure  13b  and  13c.  Hence, 
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It  represents  a  redistribution  of  the  energy  from  uf2  to  v'2  and  w*2 .  In  fact,  for  v*2  and 

w*2t  it  is  the  only  major  energy  input.  For  u'2,  it  balances  with  the  negative  production 

term,  whereas  for  u/2,  the  balance  is  with  viscous  dissipation  and  the  total  convection. 

It  is  instructive  to  compare  the  energy  redistribution  obtained  by  the  present  LES  and 
those  by  Reynolds  Averaged  Navier-Stokes  (RANS)  simulations.  In  the  RANS  simula¬ 
tions  of  Im  et  ai  (2002),  three  Reynolds  stress  models,  GL,  GL-CL  and  SSG  model,  along 
with  the  conventional  &  —  t  model,  are  used  to  compute  impinging  and  counter-current 
jets.  It  is  shown  that  in  the  strain-dominated  flows,  the  overall  predicative  capability  of 
the  Reynolds  stress  models  depends  crucially  on  the  correct  modeling  of  the  pressure 
strain  correlation,  which  is  part  of  the  redistributive  velocity  press ure- gradient  term. 
Thus  in  section  5  we  will  examine  the  eddy  viscosity  and  the  pressure  strain  correlation 
from  the  present  LES  and  compare  them  with  the  corresponding  the  Reynolds  stress 
models. 

41.  Fluctuations  in  local  s  —  n  coordinates 

Using  the  local  s  —  n  coordinates,  the  profiles  of  u'2,  and  w'2  at  s/D  =  0,  0.2,  0.8t  L6 
are  plotted  in  figure  14.  The  local  boundary  layer  thickness  at  the  corresponding  locations 
are  also  marked  on  the  right*  Note  that  because  of  the  symmetry  at  the  stagnation 
streamline,  the  profile  of  ut2  is  actually  shown  in  figure  14a  as  uj2  (at  $  —  0),  whereas 
the  u'2  at  s  —  0  in  figure  14a  is  actually  v'2.  This  correspondence  holds  only  along  the 
stagnation  streamline* 

Near  stagnation  point  s  —  0,  the  magnitude  of  streamwise  fluctuation  is  small  because 
it  corresponds  to  the  velocity  fluctuation  in  the  axial  direction  of  the  strongly  stretched 
vortices.  At  locations  s  >  0,  however,  the  streamwise  vortices  become  embedded  in 
the  mean  boundary  layer  shear  flow,  therefore  capable  of  generating  large  streamwise 
momentum  exchange  in  the  wall  normal  direction.  This  is  demonstrated  by  the  increase 
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Figure  13,  Turbulence  budget  along  the  stagnation  streamline,  a)  u ,2,  — o — :  production* 

- :  dissipation t . :  mean  convection,  - ttirbulent  convection, - :  velocity 

pressure  gradient  correlation,  b)  v'2,  — o — production,  - - :  dissipation,  - - :  mean 

convection,  — - :  turbulent  convection, - :  velocity  pressure  gradient  correlation,  c)  w*2t 

. ;  dissipation, . :  mean  convection  +  turbulent  convection, - :  velocity  pressure 

gradient  correlation. 
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of  the  magnitude  of  uf  as  one  moves  downstream.  Also  notice  that  the  peak  of  u'2 

measured  in  local  boundary  layer  thickness  is  getting  closer  to  the  wall  at  downstream, 

showing  the  characteristics  of  a  transitional  boundary  layer. 

The  ut2  shows  its  highest  peak  value  at  s  =  0  where  the  vortex  stretching  is  strongest, 
and  then  decreases  downstream.  It  becomes  significantly  smaller  than  the  st teamwise 
fluctuation  u[2  at  the  downstream  locations,  indicating  that  the  turbulence  there  is  mainly 
dominated  by  strong  shear.  Moreover,  figure  14b  also  demonstrates  the  wall  blocking 
effect  extending  to  a  distance  on  the  order  of  turbulence  integral  scale,  which  explains 
the  peaks  of  u'2  outside  of  the  boundary  layer. 

For  the  spanwise  fluctuation  wf2}  figure  14c  shows  that  its  amplitude  also  decays  mono- 
tonically  downstream  from  the  stagnation  point.  As  mentioned  before  ,  the  strong  mean 
flow  straining  and  wall  blocking  effects  decrease  as  one  moves  downstream,  and  the  flow 
gradually  changes  to  a  shear-dominated  type.  In  such  cases,  the  redistribution  term  plays 
a  much  smaller  role  in  transferring  the  turbulence  kinetic  energy  from  uf 2  to  wa,  yet  there 
is  no  production  term  for  wf2,  The  combination  of  these  factors  determines  the  trend  of 
the  decreasing  wa  in  stream  wise  direction.  Another  feature  of  the  w/2  profile  is  that 
besides  the  main  peak  inside  the  boundary  layer,  a  second  peak  develops  outside  the 
boundary  layer  in  the  leading  edge  region  (s/D  <  0.8).  We  believe  that  this  is  yet  an¬ 
other  indication  of  the  turbulence  anisotropy  resulting  from  the  existence  of  the  strong 
streamwise  vortices.  Far  from  the  leading  edge,  the  free-stream  turbulence  is  isotropic, 
and  the  turbulence  vortices  are  orientated  randomly  with  no  preferred  direction.  As  it 
approaches  the  stagnation  point,  turbulence  anisotropy  develops  because  of  the  strong 
stream  wise  vortices  emerging  from  the  background  turbulence.  If  these  vortices  were 
perfectly  stationary  then  they  would  appear  in  the  mean- flow,  but  they  are  induced  by 
free- st.  re  am  turbulence  and  thus  form  and  decay  and  move  about  as  well.  One  of  the 
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FIGURE  14,  The  profiles  of  tangential,  normal  and  span  wise  velocity  fluctuation  at  different 
stre&mwbe  locations.  The  corresponding  laminar  boundary  layer  thickness  £0.99  is  marked  on 
the  right. - :  5  =  0,  ---  — :  s  —  0.2 D, - — :  s  =  0.8i>,  and . :  s  =  1.6. D. 
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Figure  15.  The  profiles  of  turbulent  Reynolds  stress  ti'  u*n  at  different  stream  wise  locations. 

The  corresponding  laminar  boundary  layer  thickness  60,99  is  marked  on  the  right. - s  —  0, 

- :  s  =  0.2 D,  —  -■ — :  s  =  0.8D,  and  * -  :  3  =  L6 D. 

characteristics  of  the  streamwise  vortices  is  that  the  spanwise  velocity  w  changes  sign 
across  the  vortex  section,  and  thus  its  room -mean- square  value  has  a  two-peak  structure 
with  the  lowest  point  in  between  corresponding  to  the  vortex  center.  The  two  peaks  are 
not  of  equal  strength,  since  the  magnitude  of  w*  is  increasing.  When  these  streamwise 
vortices  become  strong  and  frequent  enough,  this  particular  pattern,  super-imposed  on 
the  background  turbulence  field,  gives  rise  to  the  second  peak  in  the  overall  tu'mj . 

The  only  significant  shear  stress  in  the  present  flow  is  uV  and  it  vanishes  on  the 
stagnation  streamline  due  to  the  symmetry  in  velocity  u  and  v.  In  figure  15  the  profiles 
of  are  plotted  in  local  s  —  n  coordinates.  There  are  two  distinct  peaks  for  each 
curve,  one  being  negative  outside  the  boundary  layer,  and  the  other  positive  inside  the 
boundary  layer.  Moving  downstream  from  the  stagnation  point,  one  finds  the  amplitude 
of  the  outer  peak  decreases  whereas  the  amplitude  of  the  inner  peak  increases.  The  inner 


positive  peak  of  represents  a  counter- gradient  hence  nonlocal  transport  of  inonien- 
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turn,  and  this  is  believed  to  be  a  result  of  vortex  rebounding  at  the  wall.  Essentially, 
mean  convection  brings  the  stretched  turbulence  eddies  from  free- stream  into  the  lami¬ 
nar  boundary  layer  and  aligns  the  vortices  in  stream  wise  direction.  As  these  streamwise 
vortices  approach  the  wait,  the  no-slip  wall  boundary  condition  enforces  two  thin  regions 
of  large  vorticity  (with  opposite  sign  to  the  approaching  vortices)  to  be  generated  be¬ 
tween  the  wall  and  the  approaching  vortices.  When  the  streamwise  vortices  are  strong 
enough,  the  thin  underneath  region  can  be  lifted  up  and  forms  a  dipole  vortex  with  the 
primary  vortex.  This  newly  formed  dipole  vortex  induces  a  positive  normal  velocity  and 
the  vortex  bounces  back  from  the  wall.  The  mechanism  of  rebounding  vortex  has  been 
previously  studied  in  detail  using  a  pair  of  counter  rotating  vortices,  or  a  dipole  vortex  * 
colliding  with  a  wall  under  self-induced  velocity  (Orlandi  1990;  Carnevale  et  ai  1997), 
Here  in  the  leading  edge  region*  the  vortices  are  convected  by  the  mean  flow  (thus  a 
vortex  pair  is  not  needed},  but  the  vortex- wall  interaction  is  similar.  In  terms  of  fluctu¬ 
ation  velocity,  when  turbulence  is  brought  into  boundary  layer  from  the  free- stream,  it 
creates  a  region  of  u*s  >  0,^  <  0  within  the  boundary  layer  which  is  responsible  for  the 
overall  negative  However,  because  of  the  vortex  rebounding,  the  normal  velocity 

u'n  changes  sign  in  the  near  wall  region  and  consequently  generates  the  positive  value  of 


4,8.  Thermal  fluxes  u'T'  and  v^T1 

The  profiles  of  the  thermal  fluxes  u'T7  and  u^T*  are  shown  in  figure  16  and  17  at  different 
streamwise  locations.  While  the  tangential  flux  u^T'  increases  as  one  moves  downstream, 
the  peak  value  of  the  positive  normal  flux  u^T'  decreases.  Near  the  stagnation  point,  the 
u'nTf  is  predominantly  positive  due  to  the  same  free-stream  turbulence  impingement  that 
gives  rise  to  the  negative  shear  stress  Different  from  that  of  the  however,  is 


the  peak  location  of  the  u*nTf .  Comparing  figure  17  and  15,  one  can  see  the  positive  peak 
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Figure  16.  The  profiles  of  tangential  heat  flux  u*tTf  at  different  streamwise  locations.  The 

corresponding  laminar  boundary  layer  thickness  is  marked  on  the  right.  - :  $  —  0, 

$  =  0.2 D, - :  $  —  0.BD,  and  9  =  L6D, 


of  u^T*  is  almost  five  times  closer  to  the  wall  than  the  negative  peak  of  We  believe 

this  difference  in  the  peak  locations  between  the  thermal  flux  and  the  shear  stress  explains 
the  greater  sensitivity  to  the  free-stream  turbulence  for  the  heat  transfer  rate  than  the 
wall  friction  coefficient,  as  observed  in  figure  4.  Close  to  the  wall  at  downstream  locations, 


ufnTf  changes  sign,  indicating  the  similar  counter- gradient  feature  in  the  thermal  fluxes. 
The  mechanism  can  be  understood  in  a  similar  way  through  the  vortex  rebound  as  for 


discussed  in  the  last  subsection.  The  profiles  of  the  production  terms  for  ufTr  and 
v'Tf  are  shown  in  figures  18  and  19.  Along  the  streamwise  direction,  the  u*T*  increases 
but  v'Tr  decreases.  The  decaying  of  the  production  in  vfT*  may  be  attributed  to  the 
difference  in  the  peak  locations  of  the  fluctuation  velocity  v*  and  temperature  Tr.  The 
former  shifts  away  from  the  wall  at  downstream  while  the  latter  remains  close  to  the  wall 
in  the  entire  streamwise  range. 
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Figure  17,  The  profiles  of  wall  normal  heat  flux  ufnTr  at  different  streamwise  locations.  The 
corresponding  laminar  boundary  layer  thickness  5q_$$  is  marked  on  the  right,  — — — :  5  =  0* 
— $  —  Q.2.D*  — *■ — :  s  —  0  BDt  and . :  $  =  1,6 D. 

5.  Eddy  viscosity  and  pressure-strain  correlation 
As  mentioned  in  the  introduction,  modelling  turbulence  effects  in  strain  dominated 
flows  remains  a  challenging  problem.  One  of  the  primary  goals  of  the  present  LES  study 
is  to  provide  useful  information  for  this  effort.  Conventional  two  equation  models  do  not 
address  the  anisotropy  of  stagnating  turbulence,  and  the  use  of  linear  constitutive  model* 
u'u'  =  -2&*rS;j  +  2/3HJJ,  results  in  severe  over  prediction  of  the  turbulence  production. 
In  order  to  recover  the  correct  turbulence  production*  one  remedy  (Durbin  1996}  based  on 
realizability  consideration  introduces  an  upper  bound  on  the  eddy  viscosity*  vr  ^  ji\s\ 1 
where  |5|  is  the  magnitude  of  the  strain  rate.  In  practice,  the  constraint  can  also  be 
expressed  in  terms  of  turbulence  time  scale  T  by  using  formula  ur  —  C^kT.  For  example* 
T  —  min{“,  ^ } *  where  a  =  0.6,  has  been  incorporated  into  standard  k  —  c  model 

to  compute  turbine  blade  heat  transfer  (Medic  &  Durbin  2002).  For  the  realizability 
condition*  figure  20  shows  the  normalized  turbulence  eddy  viscosity  vj  (not  the  SGS 
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Vu'T1  x  10a 

FIGURE  IS.  Profiles  of  turbulence  production  for  u'T*  along  wall  normal  direction, 
s  —  0, - ;  s  —  0,2 D, - — :  s  =  0 -.8D,  and  •  •-  ***:  $  —  1.60. 


?W'  x  IQ3 

FIGURE  19.  Profiles  of  turbulence  production  for  i/T'  along  wall  normal  direction. - ; 

s  =  0,  — — ;  s  —  0.2D, - :  s  =  0.8£?f  and  - is—  1.6 D. 

eddy  viscosity)  and  the  corresponding  along  the  stagnation  line  based 

on  the  present  LES  results.  The  effective  vr  is  obtained  through  V  —  2 with 
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Figure  20,  Turbulence  eddy  viscosity  vt  (not  the  SGS  vr). - :  LES  results, - 

k  —  e  model,  :  realizability  constraint. 


the  turbulence  production  V  and  strain  rate  Sij  computed  from  the  LES  results.  As 
a  comparison,  the  eddy  viscosity  obtained  using  the  conventional  k  -  e  model  pjjt  — 
C^Ar2/e,C^  —  0.09  is  also  plotted  .  Away  from  the  stagnation  point,  the  value  of  t^rk  is 
about  five  times  of  ur-  Close  to  the  stagnation  point,  decreases  and  approaches  to 
a  constant  value  about  two  times  of  ut.  Because  of  this  elevated  eddy  vicosity,  the  k  -  € 
model  over  predicts  the  turbulence  production  in  stagnation  point  flows.  On  the  other 
hand,  the  eddy  viscosity  vtr  tracks  the  vr  quite  well  and  provides  a  useful  upper  bound 
for  vt,  indicating  the  effectiveness  of  the  realizability  constraints. 

In  second-order  Reynolds  stress  models,  the  pressure- strain  correlation 


p  dxj  dii 


(5.1) 


is  the  primary  means  to  redistribute  the  kinetic  energy  among  different  components  of 
the  Reynolds  stresses.  The  current  LES  results  can  be  used  to  examine  the  particular 
pressure-strain  correlations  in  different  Reynolds  stress  models.  Here  three  of  such  turbu¬ 
lence  models  are  chosen  that  were  used  in  Im  et  af.  (2002):  the  GL  (with  wall  reflection 
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Figure  21.  The  pressure-strain  correlation  (no  summation  over  a)  along  the 

stagnation  streamline,  comparison  between  the  LES  results  and  different  Reynolds  stress  models, 
—  o  — :  LES,  -  GL  model,  — - :GL-CL  model, - :  SSG  model. 
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terms),  GL-CL  and  SSG  models.  Figure  21  shows  the  comparison  of  the  along  the 
stagnation  streamline  between  the  LES  results  and  the  turbulence  models  (see  Appendix 
B  for  the  expressions  of  <fiaQ  in  these  models),  First  notice  that  in  the  region  close  to  the 
stagnation  point  the  GL  model  predicts  all  three  components  of  4>aa  with  the  wrong  sign. 
Between  the  GL-CL  and  the  SSG  models,  the  GL-CL  is  significantly  better  in  predicting 
<£ii  and  033;  not  only  the  sign  is  more  consistent,  the  magnitudes  are  also  comparable  to 
the  LES  results.  The  LES  data  show  the  02 2  is  smaller  than  the  other  two  components, 
this  is  however  not  apparent  in  the  model  predictions.  Both  GL^CL  and  SSG  models 
severely  overpredict  with  the  GL-CL  model  having  the  largest  error  albeit  a  more 
consistent  sign.  Another  feature  about  the  GL-CL  model  is  that  the  peaks  of  the  0Qa  are 
located  at  a  larger  distance  from  the  wall  than  the  LES  results.  These  findings  are  con¬ 
sistent  with  the  observations  made  by  Im  et  ai  (2002)  from  the  RANS  calculations  and 
suggest  that  the  Gl^CL  model  is  the  overall  most  effective  among  the  three  turbulence 
models  for  stagnation  point  type  or  strain  dominated  flows. 

6*  Concluding  remarks 

In  this  paper,  the  stagnation  point  flow  and  heat  transfer  at  an  elliptical  isothermal 
leading  edge  in  the  presence  of  free-stream  turbulence  is  investigated  using  LES.  Along 
the  stagnation  streamline,  the  turbulence  intensity  in  streamwise  {x}  and  spanwise  (z) 
directions  increases  whereas  the  intensity  in  transverse  (t/)  direction  decreases.  Very  close 
to  the  wall,  the  streamwise  fluctuation  is  reduced  by  the  wall  blockage  effect,  and  its  en¬ 
ergy  is  transferred  to  the  other  two  components.  As  a  result,  the  spanwise  turbulence 
intensity  becomes  the  most  dominant  among  the  three  components  in  this  region.  The 
most  important  Reynolds  stress  budget  terms  for  the  stagnating  turbulence  are  the  pro¬ 
duction  and  the  velocity  pressure  gradient  correlation.  In  particular,  the  latter  dominates 
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the  development  of  the  transverse  and  spanwise  normal  Reynolds  stresses,  and  hence  the 

turbulence  anisotropy.  The  boundary  layer  under  tree-stream  turbulence  is  found  to  be 

in  non-equilibrium  and  pre- transitional  The  normal  turbulence  heat  flux  is  shown  to 

have  a  peak  location  much  closer  to  the  wall  than  that  for  the  Reynolds  shear  stress, 

explaining  the  greater  sensitivity  to  free-stream  turbulence  for  the  heat  transfer  than  the 

momentum  transfer  in  stagnation  point  flows. 

Intense,  streamwise  vortical  structures  generated  at  the  leading  edge  due  to  strong 
vortex  stretching  are  found  to  be  the  direct  cause  of  the  heat  transfer  enhancement. 
These  vortices  have  a  spanwise  dimension  about  2^3  times  the  local  boundary  layer 
thickness  and  they  move  laterally  at  a  speed  comparable  to  the  local  spanwise  fluctuation 
velocity.  The  enhancement  of  the  heat  transfer  obtained  by  the  present  LES  agree  well 
with  the  corresponding  experimental  measurements. 

The  LES  results  are  used  to  examine  the  eddy  viscosity  obtained  from  the  conventional 
k  —  e  model  and  from  realizability  conditions.  The  k  —  c  model  is  found  to  produce 
excessive  turbulent  eddy  viscosity  in  stagnation  point  flows,  but  the  realizability  condition 
provides  an  adequate  upper  bound  for  the  eddy  viscosity.  Furthermore,  the  pressure- 
strain  correlations  obtained  by  the  present  LES  are  compared  with  those  from  three 
Reynolds  stress  models  (GL,  GL-CL  and  SSG).  Among  the  three  models,  it  is  found  that 
the  Gl^CL  model  provides  the  most  effective  pressure  strain  correlations  in  stagnation 
point,  or  strain  dominated  flows.  It  is  hoped  that  the  LES  data  will  be  useful  in  developing 
improved  models  of  pressure-strain  correlations  for  strain- dominated  flows. 
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FIGURE  22.  Comparison  of  boundary  layer  velocity  profile  at  the  leading  edge.  Re  =  IQ5, 
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Appendix  A,  Validation 

Laminar  compressible  boundary  layer  profiles  at  the  leading  edge  allow  self- similar 
analytic  solutions  when  viscosity  varies  as  a  linear  function  of  temperature  (Reshotko 
&  Beckwith  1957),  Figures  22  to  23  show  the  velocity  and  enthalpy  profiles  for  Pr  —  1 
for  flow  over  an  elliptical  leading  edge.  The  wall  temperature  is  set  to  be  twice  the  total 
temperature  of  the  incoming  flow,  i.e.  Tw/Tq  =  2.  For  different  Mach  numbers,  Ala  — 
0.1  and  A/a  -  0.8(not  shown  here),  excellent  agreement  in  the  velocity  and  enthalpy 
profiles  is  obtained  between  the  numerical  and  analytical  solution.  Further  comparisons 
at  difleren  flow  conditions  can  be  found  in  Xiong  (2004),  To  validate  the  code  in  unsteady 
computations,  we  choose  to  compute  the  boundary  layer  receptivity  to  ambient  sound  at  a 
blunt  leading  edge.  Receptivity  is  defined  as  a  process  by  which  external  flow  disturbances 
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Figure  23,  Comparison  of  boundary  layer  enthalpy  profile  at  the  leading  edg e.Re  =  10s, 
M  =  0*15,  Tw/To=2,Oi  Pr  =  1,0, - :  computation ,  — — :  similarity  solution 


Figure  24.  Comparison  of  the  wall  vorticity  iuwau- - — ,  present  computation, 

result  of  Coll  is  (1997), 

are  converted  into  instability  waves  (Morkovin  1969),  For  the  flat  plate  boundary  layer 
flow,  the  sound  receptivity  refers  to  the  generation  of  T-S  instability  waves  inside  the 
boundary  layer  by  free-strearn  acoustic  waves.  We  use  the  same  flow  configuration  as  in 
Lin  (1992)  and  Collis  (1997)  for  a  compressible  boundary  layer  on  a  flat  plate  with  a 
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FIGURE  25.  St  ream  wise  velocity  profile  at  x  —  2.783. - T  present  computation, - , 

result  of  Collis  (1997). 

super  ellipse  leading  edge.  However,  unlike  the  usual  numerical  receptivity  studies,  which 
is  based  on  linearized  governing  equations  about  a  base  flow,  we  use  the  full  nonlinear 
N-S  equations  to  compute  both  the  base  and  the  disturbed  flows.  The  evolution  of  the 
disturbance  is  obtained  by  subtracting  the  base  solution  from  the  instantaneous  solution. 
In  the  present  computation,  the  mean  flow  is  judged  to  be  steady  after  the  residue  has 
dropped  by  8  orders  of  magnitude  from  its  initial  value.  The  mean  flow  wall  vorticity 
ujw  is  shown  in  figure  24  and  it  is  in  very  good  agreement  with  those  from  Collis  (1997) 
In  figure  25,  the  streamwise  velocity  profile  at  x  —  2.783  is  shown  as  a  function  of  the 
Blasius  variable  rfo,  defined  as  -  yy/Rej{x  -f  1).  Again,  the  present  result  agrees  well 
with  Collis  (1997). 

Once  the  base  flow  is  obtained,  acoustic  waves  axe  introduced  through  the  inflow 
sponge  at  a  frequency  w  “  3.31 2  and  amplitude  A  =  0.001  (see  Xiong  (2004)  for  details), 
With  the  mean  flow  Mach  number  M  —  0.1,  the  downstream  acoustic  wave  length  is 
Aa  =*  20.833, 
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Figure  26.  Disturbance  vertical  velocity  fields  of  sound  receptivity  on  a  fiat  plate  with 
a  super-ellipse  leading  edge.  Contour  levels:  Vmax  —  1-17  X  10^ 4 ,  Vmm  —  “3*90  x  10  \ 
Au  =  8.16  x  10_o* 


Figure  27.  Local  maximum  of  the  amplitude  of  the  T-S  wave  based  on  the  vertical  disturbance 

velocity,  - — - :  present  computation,  ------ :  result  of  Collis  (1997), - :  result  of  Lin 

(1992). 

The  T-S  wave  field  inside  the  boundary  layer  induced  by  the  sound  wave  is  plotted  in 
figure  26  using  the  vertical  velocity  component.  Notice  the  overall  disturbance  fiow  field 
in  the  boundary  layer  contains  not  only  the  excited  T-S  wave,  but  also  the  incoming  and 
scattered  acoustic  waves. 

After  subtracting  the  acoustic  components  from  the  total  disturbance  solution  follow- 
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mg  Wlezien  (1994  )T  the  amplitude  of  the  T-S  wave  is  shown  in  figure  27,  The  results  from 
Lin  (1992)  (for  incompressible  flow)  and  from  Collis  (1997)  are  also  shown.  Again,  we 
can  see  that  the  computations  have  followed  similar  shapes,  with  some  small  quantitative 
differences.  Given  the  great  sensitivity  of  the  growth  rate  of  the  T-S  wave,  the  difference 
is  considered  to  be  acceptable,  and  the  overall  agreement  is  satisfactory. 


Appendix  B.  Pressure  strain  correlation 


Let  <(>Qa  denote  the  pressure  strain  correlation  (no  summation  over  a)  in  the  equations 
for  u2,u2  and  w2t 


p'  9ua 
P  dxa' 


(Bl) 


and  k  and  c  denote  the  turbulence  kinetic  energy  and  dissipation.  Along  the  stagna¬ 
tion  streamline,  the  <pQQ  in  the  GL,  GL-CL  and  SSG  Reynolds  stress  models  takes  the 
following  form  : 

CL  model : 


0n  [a,(S*  -  +  2 a, .IF/,]  -  y[(2Pn  -  P22 )  -  2^(2 Pn  -  P2 i)fjp2) 

fat  =  [ait?  -  -  yI(2fta-fli)  +  &.(2Pki  -  P22)fn]  (B3) 

<fc»3  -  ~l  [ai(i?  -  /n]  +  y  [(Pm  +  P22)  -  {hw(2Pn  ~  P22)f»]  (B4) 

GL-CL  model'. 

J.  e  i  f~  2lv  ,  ^  “2 1  1  ^(2Pll  -  P22) 

011  =  [Oi(u2  -  -k)  +  2alufu2/n] - - - 

~fn [2(72 uj  +  l’2JPn  +  2y2wP22  +  2Jf(2^-^)gn)  (B  5} 

,  f  f  r 2  2 1  1  o2(2P,2-Pn) 

0 22  =  -  ^k)  -  aiwU2fn] - - - X - — 

+  /n[( l2w  +  72to)Pll  +  72^22  H - ^  (B6) 

033  =  [ai(u/z  -  | k)  -  alwu3/n]  +  Q^Pl‘3+  Pz^ 
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+/4{72i*/  +  7U^ii  +  l2wPz2  +  ^2tlJ  Sii) 


(B7) 


where 

p  -5  91/  p  ^  dv  du  dv 

pu  =  —  U2  — ,  P22  —  —  v2  -r— ,  5ii  =  -5— ,  522  —  "q—  t  In  =  - 

ax  ay  ax  oy  an  e  x* 


CBS) 


The  coefficients  are  o\  —  1.8,  <*3  =  0.6,  ctiw  —  0.6,  =  0.3,  72™  =  0.08,  Y2w  ~ 

0.1,  j2w  =  0.4,  an  =  2.5.  The  xn  is  the  normal  distance  to  the  wail. 

SSG  model: 


0n  =  -(«rie  +  *iP)&ii  F<r2€(b2n  -  |n4  +  (o3  -  <r;ni/a)A:S|i  + 

-<f4^(25n5n  -  622^22)  (B9) 

022  —  — (0TC  FajP)b22  +  az€{b22  -  -II&)  -4-  (^3  -  &ztl)I2)k$22  + 

-(74^(2622522  -  6n  5ll)  (B  10) 

033  -  “(fflt  +  crjp)633  +  036(633  —  ~nh)  —  -<74^(6n5n  -f  622522)  (Bil) 


where 


P  =  Pn  +  Pn,  n6  =  6?i-  {B 12} 

The  coefficients  are  t7]  =  3.4,  dj  =  1 .8, <J2  —  4.2,  d3  =  0.8,  d£  =  1.3,  d4  =  1.25. 
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